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Introduction 


This paper is meant as a brief survey of spectral methods, particularly 
as applied to partial differential equations. No attempt Is made to cover any 
topic In depth, but rather to present a general outline to scientists 

unacquainted with the subject. 

He begin by defining the basic spectral algorithms, emphasising 
collocation and discussing the main advantage of the method, namely the 
infinite order of accuracy attained in problems with smooth solutions. 

tn the next section ve present examples of theoretical numerical analysis 
of spectral calculations. These are elementary proofs for simple problems, 
but still may be taken as representative of more sophisticated results. 

„e conclude with an application of spectral methods to transonic flow. 
The full potential transonic equation is among the best understood among 
nonlinear equations; although there are few analytic solutions, there are many 
efficient finite difference codes. It is, moreover, of great engineering 
interest because In spite of all the simplifying assumptions introduced, it 
fits experimental data quite well. 

It was very Interesting to see what a spectral method could achieve In 
this problem, especially since there are several mathematical points not yet 
covered by any theory: nonlinearity, singularities, shock waves and entropy 
conditions. It turned out that the results are very satisfactory, and the 
algorithm as a whole is as efficient and as accurate as the best finite 
difference schemes, with the bonus of reduced memory requirements. 

Finally, a few words about the special functions mentioned In the 

paper. Only the Chebyshev polynomials 
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T (x) = cos(n cos 
n 


-1 


(x)) 


are considered in any detail, 
be taken for granted; further 
[2] and in the treatise [3] . 


Various identities and quadrature formulas will 
information can be found in references [l] and 


1. Spectral Methods - How and Why? 

. „ to find an unknown function u, satisfying some 

Assume one has to tina an 

differential equation. 

A spectra! method of solution starts * expanding u in a series oE 
eigenfunctions of a Sturm-Uouville problem. Then, using orthogonality and 
various Identities among sue, special Eunctlons one M y define approbations 
to the derivatives of u, and employ those to compute u. m practice, the 
eigenfunctions „1U usually he trigonometric functions or orthogonal 

polynomials • 

AS a simple example, consider the wave equation with periodic Initial 

data: 


= U x 


( 1 . 1 ) 


t > 0 

* 

u(x, t=0) = 4>(x); <K X ) = + 2lt) 


-i rnitrier series rather than a 
The solution will be periodic, suggesting <. 

polynomial expansion. Let 


N-l 


u N (x) = I a k (t)e 
N i n ^ 


ikx. 


( 1 . 2 ) 
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then the derivative is given by 


3u 


N-l 


(1.3) 


I Ik^Ct)- 


ikx 


3x 


k=0 


and the original problem (1.1) reduces to 


(1.4) 


da k 

dt“ = ik a k’ 


0 < k < N. 


This is now a system of uncoupled 
coefficients a k . The only input 
initial values of a k (0) . These are 
^>5 they are defined as integrals 

(1.5) V 0) = 27 


ordinary differential equations for the 
..eeded for solving (1.4) is the set of 
the Fourier coefficients of the function 

j 2 * e - lkx ♦wax, 

0 


or, in a more practical fashion, as sums 


( 1 . 6 ) 


N-l -ik 


2i\l 


Tc 


(0) = i I e 


1=0 



which can be computed in 0(N log N) operations by means of the fast Fourier 
transform (FFT)* 

Ar. alternative approach respiting in (1.0, hut generalitable to 
arbitrary orthogonal series is detUn's method: Substitute (1.2) into 

0.1), nuitlply (1.1) by for k - 0,1 »-l and integrate over the 

period. (For another system of eigenfunctions, substitute an approximating 
sum, multiply by the eigenfunction and the corresponding weight function, and 
integrate over the interval of orthogonality.) 


- 4 - 


We now impose boundary condition on the equation: 


( 1 . 7 ) 


u t = u x> |x| <1, t > 0 

u(x, t=0) = (|»(x) 
u(x=l, t) = l(/(t). 


A Fourier series is no longer appropriate, so we expand in terms of Chebyshev 
polynomials : 


0 . 8 ) 


N 

“h ‘ J 0 \ (l;)T k (x) - 


The formula for the derivatives is more complicated 


( 1 . 9 ) 


9x 


N 


l 

k=0 


b k(t)T k (x) 


C k b k 


N 

= l 

p=k+l 
p+k odd 


P a 



( 2 

if 

k = 

0 

c v = 

1 




iC 

( 1 

if 

k > 

0 


but the final equation for the coefficients is similar to (1.4) 


- 


” 

a o 


a 0 

a l 


a l 

• 

• 

= A 

• 

• 


* 

a N 


a N 

- 


to 




( 1 . 10 ) 


d 

dt 


„ lt h . certain matrix A. The Initial walues for the coefficients can again 
be taken to be the Chebyshev coefficients of the Initial data which can be 
efficiently computed by frt. Thete la, however, an additional elation to be 

satisfied : 


( 1 . 11 ) 


l a. (t)T k (l) - *(0, 
k=0 


which represents the boundary condition. Note that (1.10) and (1 
overspecify the H unknown coefficients; this would be no Problem if infinite 

series were used, but for a calculation with a finite number of .odes some 

, . . _ satisfying only N of the N+l equations 

compromise has to be made, e.g., satisry g 

. , t , ,, u) Reserving some of the coefficients to satisfy 

(1.10), together with (1.1 U. Keservxi 

boundary condlrlon. - as done here - Is called a Lanosos fia-tM. 

Finally, consider a wave equation with variable coefficients: 


( 1 . 12 ) 


u t = (c(x)u) x , 
u(x, t=0) 
u(x=l,t) = Kt). 


< 1, t > o 


a(x) > 00 > 0. 


one should again approximate » by an Nth degree polynomial, hut the 

coefficients of c(x)u cannot usually be defined in terms of the coefficients 

of u. This forces us to adopt a different approach: 

Take Nil points * 0 , x, x„ In (-Ml- These define a unique 

polynomial of degree N which is Identical with u at the points - the 

Intetpolant of u. We now replace c(x)u by the interpolant of c(x)u - which 

since all is needed are the values of u(xj) “ and 


is readily available, 


compute its derivative to advance (1.12) in time. The boundary condition is 

satisfied by having x Q = 1 and by setting u(x Q , t) = 'J'(t). 

This procedure - which can, obviously, be applied to nonlinear equations 
too, is called a collocation method or a pseudospectral method. It is more 
general than the Galerkin and tau methods mentioned above, and boils down to 
defining the values of f'(xj) given f( Xj ), accurately for all polynomials or 
trigonometric polynomials of degree < N. This can be of course done for any 
set of points, and the corresponding operator is represented by a matrix 

multiplication 

(1.13) = l V f(X k )> 

tv 

2 

However, this is an inefficient numerical procedure, needing 0(N“) 
operations. For special sets of collocation points the matrix multiplication 
can be done by FFT in 0(N log N) operations. This is, in fact, one of the 
reasons why trigonometric functions and Chebyshev polynomials are usually 
employed in spectral calculations. Note that once the algorithm (1.13) is 
available, the expansion coefficients are no longer needed - in contrast to 

formulas (1.4), (1.10), (1*11)« 

In conclusion, we have introduced three kinds of spectral methods: 
Galerkin, tau, and collocation. We have singled out Chebyshev and Fourier 
collocation methods as the most useful for two reasons: 

a) They may be used in variable coefficient and nonlinear problems; 

b) They allow the fast Fourier transform. 


We now address the second question in the title of this section: Why use 
spectral methods? To answer this, notice that the spectral approximation 
obtained in (1.2) - (1.4) differs from the exact solution of problem (1.1) by 
the quantity 


( 1 . 14 ) 


00 


I 

k=N 


ik(x+t) 

e 


where a k are the Fourier coefficients of $. The error is certainly 
majorized by the sum of the absolute values of a^. Now, if one assumes a 
smooth 4>, i.e., one possessing continuous derivatives of all orders, it turns 
out that the coefficients decay faster than any power of k: 


( 1 . 15 ) 


- 4 ! 2 ' e- ikX *(K)dx = ^ r e_ikX ♦'<*>** 


k 2ir 


2v 


2ir 


2tt ( ik) 0 


/ e -1KX <{> (x)dx 


M 


•2tt( ik) 0 


f -ikx , (M) / . . 

/ e y (x)dx, 


(simple integration by parts; M is arbitrary). Therefore, a spectral method 
using N modes and applied to smooth functions will admit an error estimate 
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(1.16) 


| error | 


C(M) 

.,M 


> 


for any M, with some constant C(M). (In contrast, a finite difference 

c 

method with N grid points will usually have an error of the form — where 
p is some small integer such as p = 2 for a second order method.) Formula 
(1.16) is usually referred to as the "infinite order of accuracy of a spectral 
method," or, for short, "spectral accuracy." It also holds for general 
eigenfunction expansions; provided one deals with smooth functions. 

Computationally, (1.16) means one can obtain very accurate numerical 

solutions, using relatively few data points. 

We should also say a few words about nonsraooth functions. This is not an 

academic question, since most nonlinear hyperbolic systems admit discontinuous 

solutions, and even allow discontinuities to evolve in time from smooth 

initial data. For example, consider the Euler equations of gasdynamics, 

which produce solutions with shock waves, contact discontinuities, and 

rarefaction waves. In this case one cannot expect high accuracy - rather, the 

Gibbs phenomenon which occurs at discontinuities will produce an oscillating 

error component which does not vanish as N ». However, proper treatment of 

this problem may filter out the noise and produce good approximations "away 

from shocks" [4], In fact, the Gibbs phenomenon itself may serve as a shock 

accurately pinpointing sharp transitions in the solution. 


locator, 


Although there exist several satisfactory spectral calculations of nonsmooth 
solutions [5, 6, 7], this area Is very much In need of a firm theoretical 


basis. 


2. Proofs 

In this section we present two examples of analysis of spectral 
methods. The problems we treat are of the form 


and are replaced by numerical approximations 

3u N 

atT T! V 

We present elementary convergence proofs, showing that the numerical results 

actually approach the unknown functions sought. 

We use the Lax equivalence theorem, which states that a scheme which is 

consistent and stable Is convergent. Consistency Is the following 


property 
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It expresses the fact that the spectral operator approaches the differential 
operator as the number of modes increases, and will be taken for granted, 
under the assumption of spectral accuracy. What we prove is stability, i.e., 
estimates of the form 


Ju N (x,t)n < c nu N (x, t=o)n, o < t < t, 

which should hold in some norm, with constants C which may depend on T, but 
not on N, the number of modes employed. 


2.1 Fourier Collocation 

We solve the wave equation with periodic boundary conditions 


u t = V 


( 2 . 1 ) 


t > 0 

u(x, t-0) = $(x); $(x + 2 tt) = 4><x), 

2*j 


by Fourier collocation, using the points x^ , 0 < j < N; we assume that 

N is even, N = 2M. The approximate solution u is a trigonometric 

polynomial of degree N, representable as a linear combination of the 


functions: 


( 2 . 2 ) 


1, cos(x), cos(2x) , . . . ,cos(Mx) 

sin(x), sin(2x) , ••• ,sin((M-l)x) 
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which satisfies: 


(2.3) 


3u 3u 

3tT (x j J = W (x j ) = l D jk “n^ X j ^ ' 


From (2.3) one can deduce the stability of this numerical scheme, namely: 


(2.4) 


2 2 
0<j<N " J 0<j<N J 


One way of proving this is by explicitly computing the collocation derivative 
matrix Dj^, which turns out to be antisymmetric. Therefore, when (2.3) is 
multiplied by u^(xj,t) and summed, one obtains 


(2.5) 


It J l U N (x i’ c) = £ D ik VV t)u N (x k ,t:) = 0 

0< j<N 3 0<j,k<N 3 N J N * 


which implies (2.4). 

The second method to establish stability is based on a different 
interpretation of the sums in (2.5): 


( 2 . 6 ) 


9u v 


3 1 r 2, * r , . N , N 

3t 7 4 u N (x i»t) = l u (x t;^— (x t) ( 

dt 1 0<j<N N 3 0<j<N N 3 3x 3 


An examination of the functions in (2.2) shows that one may replace the second 

^ U N 

sum by an integral, because for any combination appearing in u. T - — the 

N oX 

trapezoidal rule with N points is exact: 


(2.7) 


2lT „ 

/ f (x)dx = l f(x ) , 

0 0<j<N J 




Thus, 


gfifi 


(2.8) I V X j’ C) 

dt Z 0< j<N J 


_N f 2lt u ( x ) — - dx = 0 (by periodicity) 
2u q N °x 


leading again to (2.4), 


2 . 2 Chebgshev Collocat ion for an Initlal-Bourn^a^^ 
The differential equation to be solved is: 


(2.9) 


u = u 

t X 


< 1, t > 0 


u(x, t=0) = <|>(x) 
u(x=l , t) = 0. 


Consider the points: 


( 2 . 10 ) 


jx 

X = cos , 


0 < j < N. 


The identity 


- 1 > j nm 

T (x.) = cos(m cos (x.)J = cos — jj- 
mi J 


show, clearly how one can use FFT <a cosine transfer™) to fit Interpolating 
polynomials to data at these points. To solve (2.9) we begin by satisfying 

the boundary condition by taking: 


( 2 . 10 ) 


u N (x 0 ) = °’ 
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while the collocation procedure is: 


( 2 . 11 ) 


3u M 3u m 


0 < j < N. 


(Note that the point x N - -1 is not a 
condition, one should set 


collocation point.) For the initial 


( 2 . 12 ) 


u N (Xj , t=0) = <f>(Xj), 


0 < j < N. 


Conditions (2.10) - (2.12) uniquely define a polynomial u N (x.t), of degree 

N-l in x, with coefficients depending on t. 

A little thought shows that (2.11) may be replaced by 


(2.13) 


3u„ 3u M 


nhero p N U> is the unique polynomial of degree K-l »hith vanishes at . 

0<j <N and takes the value 1 at x Q . In fact, since T N (x) 

can explicitly identify p N 


extrema in [-1,1] at the points x j , 


one 


(2.14) 


T'(x) 

V x> ■ “ 7 " • 


Thus one may extend formula (2.11) - equality at certain points - to the 
formula (2.13) which holds everywhere. This is done at the expense of an 
nahno^ function r(t); we stress again that P»00 is an explicitly hnown 

function. 

Multiply now (2.13) by u n 

J 1 _ V 


and integrate: 


pr: y 


1 > \ ■* " ' " * ' " ■ . - V "" ” v • K 


(2.15) 


d r l i 2 1 + x dx = f 1 u !!« 1 + _* 

dt 2 N j 2 -1 N 3X V 1 - x 

/ 1 - x 


1 + X 


V^7#i dx + ' <t) ^- p »7^ 


The first integral on the right is negative, as integration by parts shows: 


(2.16) 


2 2 

9u N 1 + x U N 1 + x 1 _ f 1 ±_ /HE 

^ l «.l«-y=y-^7==f il 2d * /l '* 

1 / 1 - X / 1 - X x=-l 


..| 2(h *r 1/2 (i - >•)■ 


The other two integrals may be replaced by sums over Xy Indeed, the 
following formula holds whenever g is a polynomial of degree < 2N 


(2.17) f 1 g — — dx - I c g(x ), (<=0 " C N “ 1 1 C 1 " *> 

-‘/IT? i-° 


Then the multiplier of T 


(2.18) 


/ u n P N 7 1 - + -% " J - C J U N (x j )p N (x j )(1 + V 
-1 / 1 - x J 


is seen to vanish, since 
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1 + x = 0 at x = X N = -1 
u N = 0 at x - x 0 = 1 

p N = 0 at the interior points Xj , 0 < j < N 


In conclusion, we find 


( 2 . 19 , /,*>♦* 


1 2/ , , „n“ 1/2 (! _ x )" 3/2 dx < 0 


and have proved stability in the form 

nu N (x,t)i < nu N (x, t=o) ii 

where the norm is defined by 

(2.20) «^» 2 " / /f?l u2(x)dX * 

A few remarks are in order about the above results. We note that special 
properties of the expansion functions are extensively used throughout. 
Orthogonality is needed, especially in the subtler form of Gauss quadrature 
formulas; for instance, formula (2.17) uses N+l points, but is accurate for 
polynomials of degree 2N-1. None of the methods used in this section would 
work for arbitrary collocation points; in particular, polynomial collocation 
at evenly spaced points has very different (and numerically undesirable) 
properties. This is, in fact, the basic reason why eigenfunctions of Sturm- 

used in spectral approximation; again, the Chebyshev 


Liouville problems are 



polynomials stand out as readily available functions, allowing general 
boundary conditions, and admitting the Fast Fourier Transform. 


3. Applications to Transonic Flow 
3.1 The Problem 

Consider the steady two-dimensional flow of a compressible, inviscid gas 
past an airfoil, with uniform conditions at infinity. Under these 
assumptions, one may take the velocity vector to be the gradient of a velocity 
potential $: 


(3.1) « = * x , v = 

The density p is determined by Bernoulli's law 


(3.2) 




+ v 


1 1 



where y is the specific heat ratio and the free stream Mach number. 

The equation one must solve, for the scalar unknown 0, expresses mass 
conservation: 


(3.3) 


(p<& ) + (p$ ) - 0* 

x x y y 


The boundary conditions for (3.3) are as follows: 




a) At infinity, the potential must satisfy 


(3.4) 


$ ~ x < > u ~ 1 


representing uniform flow with normalized horizontal speed 1. 
b) Let the wing be located on the x axis, with a shape given by 


(3.5) 


±t( 1 - x ), 


Ixl < 1, 


- a parabolic wing of (small) thickness ratio = t. We impose 
the boundary condition 


(3.6) 


9$ _ dy 3$. 
Ity ~ dx 3x 


at y = 0, |x| < 1, 


which 


approximates , to the first order In r, the exact condition 


(3.7) 


|t- 0 on 
3n 


the body y = ±t( 1 - x ) 


Because of symmetry, 


(3.8) 


li= 0 


3y 


on the rest of the x axis. (The symmetry about the x axis means, 
of course, that there is no lift force on the wing; this assumption is 
made only for mathematical convenience, i.e., solving for y > 0 


The outstanding property o£ equation (3.3) is that it changes type fro. 
eliiptic to hyperbolic as the local speed goes fro. subsonic to supersonic, 
we shall describe briefly the qualitative behavior of the flows it represents. 

is long as the Mach number M„ is small enough, the flow is smooth and 
symmetrical under the transformation x — > -x; the singularities at x i 
do not propagate into the field. The speed is either subsonic everywhere, or 
a small supersonic pocket may develop over the wing. A, H. increases over a 
critical value, a shock wave develops near the trailing edge, across which the 
speed reduces abruptly from supersonic to subsonic. Thus the solution is no 

longer symmetric , although there is symmetry in the differential equation and 

. , ..--of-rlcal flow would contain an unphysical 

in the boundary conditions. A symmetrical now 

rarefaction shock near the leading edge. 

This is clearly a case of nonuniqueness of weak solutions, to be resolved 

by an entropy eondltion. It means, from the computational point of view, that 

some entropy inequality - or equivalently, some desymmetrising prooedure - 

must be enforced, in addition to equations (3.3), (3.4), (3.6), (3.8). 


3.2 The Numerical Algorithm 

As a first step towards the solution, we notice that the discontinuity of 
H at x = ±1 requires a mesh refinement there. Since we also need buffer 
Ls i„ front of the airfoil and behind it, we discrete the problem by 
using three Chebyshev domains, as shown: 





Three Chebyshev Meshes; BC Represents the Airfoil. 


We impose continuity of $ and on the interfaces BG and CF« The various 
derivatives appearing in (3.3) are computed by collocation. After these 
obvious steps are taken, one is left with two problems: 

a) to devise a desymmetrizing algorithm; 

b) to find an effective iteration scheme for the solution of the 

(nonlinear) discrete version of (3.3). 

Streett [8] has developed an algorithm to deal with these questions, which is 
moreover applicable to lifting flow, including the exact flow tangency 
condition (3.7). We summarize now Streett's method, as applied to the present 
problem. 







In finite difference codes, desymmetrizlng is done by computing 
derivatives by upwind differencing in supersonic regions [91 or by biasing the 
density towards upwind values [101. A. spectral derivative calculation is 
nonlocal, it allows no procedure similar to upwind differencing; it is much 
more convenient to use the modified density approach. A modified density P 
is computed by the formulas: 


= p - y 6 p 


(first order) 


p = p - y(6p - e6E -1 p) (second order), 


with switches y and e defined by: 


= max 


(o, i - 4) 


z = max 


( 0 , 1 - 1 ^) 


difference operators E, 6 are as follows. 


(6f)(x 1 ) = fCXj,) “ 


(6°f)(x i ) = \ (f(* 1+1 > - f ^ x i-i^ 
(E _1 f)(x i ) - f ( x i-l.) 


where x^ < x ± < * 1+1 are three neighboring mesh points. M is the local 
Mach number and K of the order of the density jump across the shock. After 
o is computed, the equation (3.3) is replaced by 
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(3.9) 


(p$ ) + (p$ ) = 0. 

v x v x ^ y'y 


Finally, we address the iteration method. One may regard $(x,y) as the 
time independent solution of an equation 


(3.10) 


A = (p$ ) + (p$ ) = N (<&) 

3t x x y y sp 


where A. is a linear operator and N stands for the nonlinear part, 

sp 

computed spectrally. Indeed, if (3.10) possesses steady state solutions, 
these will satisfy (3.9). The operator A, which seems arbitrary at this 
point, must, in fact, be chosen very carefully to ensure convergence. The 
usual choice for the transonic equations is 


(3.11) 


a$ + = N ($>) 

t xt sp 


with variable coefficients a and £: a is nonzero at subsonic points, and 
& is nonzero at supersonic points. 

The time introduced in equations (3.10) or (3.11) is purely artificial; 
we may replace the variable t with an iteration index n. Motivated by 
(11), the following iteration scheme is produced: 


(3.12) a($ n+1 - $ n ) + e($ n+1 - $ n ) » N ($ n+1 - $ n ) + N ($ n ). 

N v X x sp sp 


The two formulas (3.11) and (3.12) are not exactly equivalent, but it is seen 

that if - $> n -»• 0 then 13 ($ n ) + 0, so we are still solving (3.9). 

sp 




it involves nonlinear 


As it stands, (3.12) is still hard to solve, as 
operators and the full matrices of spectral derivation. We overcome this 
difficulty by the method of approximate inverses, namely by replacing N sp by 
a linear operator L, which is easy to invert and is near - in some sense - 

to N sp . 

We shall take L to be the finite difference representation of 


on the spectral Chebyshev grid, with P considered as a given coefficient. 
This is a five-point operator, since five points are sufficient to define the 
derivatives involved; it is also clear that L$ - * 

From the resulting iteration scheme 


(3.13) 


a($ 


n+1 _ $ n ) + g($ n+1 - $”) = U $ n+1 " $ n ) + V* n >« 


one can readily compute * n+1 - by inverting the operator a + P - L. 

This is efficiently done by a dimensional split - two tridiagonal matrices for 
separate x and y relaxation - or by an approximate LU-f actorization of the 
five-point operator. (For additional material on iteration schemes for 

spectral operators, see (11] •) 

Using the apparatus mentioned above, the transonic flow was computed for 
the simple case of the symmetric parabolic airfoil, as well as for more 
meaningful lifting airfoil shapes [8]. The results were compared with state 
of the art finite difference schemes for the same problems, and found to 
produce the same resolution with significantly smaller number of grid 
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points. Shocks were accurately captured. Although a spectral iteration is 
more complex than the corresponding finite difference one, and usually one 
needs more iterations to converge in a spectral code, the running times for 
the two methods were comparable, and in certain cases the spectral code ran 
faster then the finite difference code. This clearly shows the advantage of 
smaller computational meshes, which are admissible because of spectral 
accuracy. 
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